Content:

Ⅰ. Introduction

The m6AConquer database offers curated quantitative m6A profiling sequencing data sets. This page provides a brief overview of the processing workflows (Figure 1) utilized to convert raw reads into the data available in this database.

Ⅱ. Data Processing

1. Data gathering

m6AConquer summarises sequencing datasets from 10 different m6A profiling techniques. Both wild type and perturbation samples are included in the database. The perturbation types mainly include the manipulation of m6A-related genes, such as the over-expression of FTO and the knockout of Mettl3. IVT control samples are also included for quantification calibration.

All the sequencing datasets are downloaded from three central database (GEO, GSA, and EMBL-EBI).

2. Quality Control

FastQ read files for NGS techniques are retrieved and downloaded accordingly. Most reads undergo initial trimming using tools like Cutadapt, Trim Galore, and Flexbar to eliminate adaptors and low-quality reads. Some reads are deduplicated through Samtools, Fastx Collapser, and UMI-tools to remove PCR duplications. Fastp is used for both trimming and deduplication.

As for Oxford Nanopore sequencing data, the raw Fast5 files are obtained from the Singapore Nanopore-Expression Project. Quality control is conducted during the basecalling process with Guppy.

The table below summarize the software and tools used in data cleaning and read mapping:

3. Alignment

For NGS sequencing data, 4 aligners (STAR, Bowtie, HISAT2, and HISAT-3N) are utilized for different techniques based on their original pipelines. Bowtie is employed explicitly for transcriptome alignment, while the other three aligners are utilized for genome alignment. The genome index of HISAT2 can be obtained from the provided link: https://genome-idx.s3.amazonaws.com/hisat/grch38_genome.tar.gz. Meanwhile, the genome indexes for STAR and HISAT-3N are constructed as per the specific requirements. All sequencing reads are aligned to the GRCh38 or hg38 human genome, and various versions of gene annotation from RefSeq, Genecode, and Ensemble are employed under their original pipelines. The codes for index building can be accessed on GitHub.

In the context of Oxford Nanopore sequencing data, the aligner minimap2 is utilized for transcriptome alignment. The Ensemble transcript reference of GRCh38 is indexed using the default parameters of minimap2.

4. Feature counting

The m6AConquer database summarizes quantitative m6A profiling results in the form of SummarizedExperiment (SE) objects. These SE objects summarize the total and m6A read coverages at a consistent range of genomic locations and gene expression levels. The consistent genomic locations consist of 2526220 adenosines (A) at the center of exonic DRACH motifs and 262485 GLORI m6A sites discovered at other genomic locations reported in the supplementary data of GLORI. The m6AConquer database also includes gene expression matrices in SE objects on 58560 genes (Transctipt-Omics).

5. Omics data frameworks

The SummarizedExperiment objects consist of four primary components: row ranges, column data, assays matrices, and metadata.

In SE objects for m6A-omics (Figure 2),

In SE objects for transcript-omics (Figure 3),

Apart from these single-omics data in SE format, multi-omics data in the form of MultiAssayExperiment (MAE) for each profiling technique are also provided. The MAE files are the integration of m6A-Omics and Transcrit-Omics SE objects. The MAE contains distinct components such as experiments, sample maps, and column data for MAE.

In MAE objects for multi-omits (Figure 4),

All three types of files (multi-omics, m6A-omics, and transcript-omics) can all be downloaded from the Download page.

5. IVT calibration

In the m6AConquer, we have incorporated in vitro transcription (IVT) control samples from GLORI, eTAM-seq, and MeRIP-seq datasets. The IVT control sequencing data from GLORI and eTAM-seq employ a beta-binomial mixture model to call m6A sites separately. Sites with a false discovery rate (FDR) below 0.05 are considered false positive m6A sites for the respective GLORI or eTAM-seq datasets. In the case of IVT control sequencing data from MeRIP-seq, m6ACali is used to identify false positive m6A sites. For other sequencing data within these three datasets, the m6A counts and total counts at their identified false positive m6A sites are set to 0.

III. Site-calling and Normalization: Beta-binomial mixture model

Beta-binomial mixture model

The two-component beta-binomial mixture model is a statistical model that assumes the data is generated from a mixture of two beta distributions. Each beta distribution represents a different underlying mechanism or subgroup within the population. The probability of observing success for each individual is assumed to come from one of the two beta distributions, and the choice of distribution is determined by a latent (unobserved) variable that follows a binomial distribution:

where $\pi$ and $1-\pi$ are the proportion of two beta distribution components, and $\alpha$ and $\beta$ are beta distribution parameters.

Expectation-Maximization (EM) algorithm

In m6AConquer, we model the m6A methylation ratios detected in each sequencing sample as a two-component beta-binomial mixture distribution representing two possible methylation states: methylated state (foreground, “$\text {fg}$”) and unmethylated (background, “$\text {bg}$”) state. The fitting of the beta-binomial mixture model follows the Expectation-Maximization (EM) algorithm. The EM algorithm is an iterative method that starts with initial estimates for the parameters and then changes between the E-step (Expectation) and the M-step (Maximization). In the E-step, the log-likelihood is computed the based on the current parameter estimates. In the M-step, the log-likelihood is maximized to update the parameter estimates.

1. Initialization

For each site with m6A read coverages, let $m_i$ and $n_i$ be the m6A and total read coverages at site $i$. Let $\{x_i = \frac{m_i}{n_i}\},\ i=1,2\ldots,N$, denotes the m6A methylation ratios at site $i$, where $N$ is the number of sites with m6A read coverage.

To obtain a closer initial estimate for $\pi$, a two-component binomial mixture distribution is fitted for $x_i$:

The fitting of this primitive binomial mixture distribution is also implemented by EM algorithm.

To obtain initial estimates for $\alpha_{\text {fg}}$, $\beta_{\text {fg}}$, $\alpha_{\text {bg}}$ , and $\beta_{\text {bg}}$, weighted method of moments is firstly used for estimation (18-21), following the weighted maximum likelihood estimation (8-17) to obtain a closer estimation. The details of these two estimation methods are introduced in the M-step.

2. E-step

In the E-step, the algorithm calculates the value of the likelihood function concerning the conditional distribution of the latent variable given the observed data and current estimate of the parameters. This is where posterior probabilities are computed:

where ${r_i}_{\text {,fg}}$ and ${r_i}_{\text {,bg}}$ denotes the posterior probability of $x_i$ belonging to the foreground and background component respectively, $m_i$ and $n_i$ be the m6A and total read coverages, (4) denotes the probability mass function of beta distribution component, and $\text B$ denotes beta function.

3. M-step

In the M-step, the algorithm finds the parameters that maximize the likelihood. By default, weighted maximum likelihood estimation (WMLE) is used to obtained the estimated parameters. If the WMLE fails to converge, weighted method of moments is used in current iteration as the estimation method.

• Update component proportion

  The update formula for component proportions are:

  $$\begin{equation} \pi = \frac{\sum_{i = 1}^{N}{r_i}_{\text {,fg}}}{N} \end{equation}$$

  $$\begin{equation} \pi_j= \begin{cases} \pi & j = \text {fg}\\1-\pi & j=\text {bg} \end{cases}\end{equation}$$

  where $N$ is the number of sites with m6A read coverage.

• Update beta distribution parameters
  • Weighted maximum likelihood estimation (WMLE)

    To avoid underflow, the log-likelihood of beta-binomial mixture distribution is calculated:

    $$\begin{equation} \begin{aligned} l &= \sum_{i=1}^{n}\sum_{j\in\{\text {fg,bg}\}}[\log \pi_j +\log {p(m_i|n_i,\alpha_j, \beta_j)}] \end{aligned} \end{equation}$$

    To find the estimates of $\alpha$ and $\beta$ that maximize the log-likelihood function (8), the Newton-Raphson method is used. The Newton-Raphson method is an iterative root-finding algorithm that uses the first and second derivatives of a function to find the values of the parameters that maximize the function. The update formula in the Newton-Raphson method are:

    $$\begin{equation} parm_{\text {new}} = parm_{\text {old}} - Hessian_{\text {inv}} *f = \begin{bmatrix} \alpha_{\text {new}}\\\\ \beta_{\text {new}} \end{bmatrix} \end{equation}$$

    $$\begin{equation} Hessian = \begin{bmatrix} \frac{{\partial^2 l}}{{\partial \alpha^2}} & \frac{{\partial^2 l}}{{\partial \alpha \partial \beta}} \\\\ \frac{{\partial^2 l}}{{\partial \beta \partial \alpha}} & \frac{{\partial^2 l}}{{\partial \beta^2}} \end{bmatrix} \end{equation}$$

    $$\begin{equation} f = \begin{bmatrix} \frac{{\partial l}}{{\partial \alpha}} \\\\ \frac{{\partial l}}{{\partial \beta}} \end{bmatrix} \end{equation}$$

    $$\begin{equation} parm_{\text {old}} = \begin{bmatrix} \alpha_{\text {old}}\\\\ \beta_{\text {old}} \end{bmatrix} \end{equation}$$

    where $Hessian_{\text {inv}}$ is the inverse of Hessian matrix (10), $\frac{{\partial l}}{{\partial \alpha}}$ and $\frac{{\partial l}}{{\partial \beta}}$ are the first partial derivatives of (6) with respect to $\alpha_{fg}$ and $\beta_{fg}$ or $\alpha_{bg}$ and $\beta_{bg}$, and $\frac{{\partial^2 l}}{{\partial \alpha^2}}$, $\frac{{\partial^2 l}}{{\partial \alpha \partial \beta}}$ and $\frac{{\partial^2 l}}{{\partial \beta^2}}$ are the second partial derivatives. $\alpha_{\text {old}}$ and $\beta_{\text {old}}$ are $\alpha_{\text {fg}}$ and $\beta_{\text {fg}}$ or $\alpha_{\text {bg}}$ and $\beta_{\text {bg}}$ of the previous iteration, and $\alpha_{\text {new}}$ and $\beta_{\text {new}}$ are the updated parameters in current iteration.

    The derivatives are calculates as follow:

    $$\begin{equation}\begin{aligned} \frac{{\partial l}}{{\partial \alpha}} = \sum_{i=1}^{N} &r_{i}(\psi(\alpha + \beta) - \psi(\alpha) \\&+ \psi(\alpha + m_{i}) - \psi(\alpha + \beta + n_{i})) \end{aligned}\end{equation}$$

    $$\begin{equation}\begin{aligned} \frac{{\partial l}}{{\partial \beta}} = \sum_{i=1}^{N} &r_{i}(\psi(\alpha + \beta) - \psi(b) \\&+ \psi(\beta + m_{i} - x_{i}) - \psi(\alpha + \beta + n_{i})) \end{aligned}\end{equation}$$

    $$\begin{equation}\begin{aligned} \frac{{\partial^2 l}}{{\partial \alpha^2}} = \sum_{i=1}^{N} &r_{i}(\psi_{1}(\alpha + \beta) - \psi_{1}(\alpha) \\&+ \psi_{1}(\alpha + m_{i}) - \psi_{1}(\alpha + \beta + n_{i})) \end{aligned}\end{equation}$$

    $$\begin{equation}\begin{aligned} \frac{{\partial^2 l}}{{\partial \beta^2}} = \sum_{i=1}^{N} &r_{i}(\psi_{1}(\alpha + \beta) - \psi_{1}(\beta) \\&+ \psi_{1}(\beta + n_{i} - m_{i}) - \psi_{1}(\alpha + \beta + n_{i})) \end{aligned}\end{equation}$$

    $$\begin{equation} \frac{{\partial^2 l}}{{\partial \alpha \partial \beta}} = \sum_{i=1}^{N} r_{i}(\psi_{1}(\alpha + \beta) - \psi_{1}(\alpha + \beta + n_{i})) \end{equation}$$

    where $\psi$ represents the digamma function, $\psi_1$ represents the trigamma function, $r_i$ are the responsibilities calculated in the E-step, $m_i$ are the m6A read coverage at site $i$, and $n_i$ are the total read coverage at site $i$.

    The parameters for foreground and background components are estimated separately.

  • Weighted method of moments (WMM)

    If WMLE fails to converge in a certain iteration, the estimated $\alpha$ and $\beta$ is substituted by those estimated by weighted method of moments (WMM):

    $$\begin{equation} \alpha_{\text {mom}} = \hat{\mu} \cdot \left( \frac{\hat{\mu}\cdot (1 - \hat{\mu})}{\hat{\sigma}^2} - 1 \right) \end{equation}$$

    $$\begin{equation} \beta_{\text {mom}} = (1 - \hat{\mu})\cdot \left( \frac{\hat{\mu}\cdot (1 - \hat{\mu})}{\hat{\sigma}^2} - 1 \right) \end{equation}$$

    where $\hat{\mu}$ is the weighted mean and $\hat{\sigma}^2$ is weighted variance:

    $$\begin{equation} \hat{\mu} = \frac{\sum_{i=1}^{N} r_i \cdot x_i}{\sum_{i=1}^{N} r_i} \end{equation}$$

    $$\begin{equation} \hat{\sigma}^2 = \frac{\sum_{i=1}^{N} r_i \cdot (x_i - μ)^2}{\sum_{i=1}^{N} r_i} \end{equation}$$

    Here, the weight $r_i$ are the responsibilities calculated in the E-step.

    The parameters for foreground and background components are estimated separately.

4. Convergence

The algorithm checks for convergence by testing if the absolute differences between the current and the previous parameter values are smaller than a certain tolerance. If not, the parameters are updated, and the next iteration starts. If yes, the process stops, and the final estimates for the parameters are returned.

Model outputs

The beta-binomial mixture model in m6AConquer outputs two statistics: m6A site probabilities and false discovery rates.

where $x$ is the m6A counts, $F(x;\alpha_{bg}, \beta_{bg})$ and $f(x)$ is the cumulative density function (CDF) and probability mass function (PMF) of the background beta-binomial component respectively, $N$ is the number of sites with m6A read coverage and $\text B$ is the beta function. The CDF of the background component is calculated through the sum of PMF, which can be implemented with a recursive algorithm.

IV. Technical Orthogonal Integration: IDR analysis

IDR analysis

To identify reliable modification sites and investigate differences between techniques, m6AConquer implements Irreproducible Discovery Rate (IDR) analysis, a method commonly used in CHIP-seq data to identify consistent peaks in replicates. In our implementation, we evaluate the consistency of identified sites between different techniques, particularly those employing orthogonal mechanisms to quantify m6A abundances. IDR analysis involves three key steps to assess the reproducibility of methylation sites in various m6A profiling techniques (Figure 5).

1. Site ranking

The overlapped sites between each pair of profiling techniques are ranked based on the value of interest. In m6AConquer, m6A site probabilities are selected as the evaluation value. Therefore, the overlapped sites are ranked based on posterior probabilities, and the ranks are used as modelling statistics in the following step.

2. Modelling

The underlying assumption in IDR analysis is that reproducibility sites should rank higher and consistently in each pair of techniques. To model the heterogeneity, a mixture framework is employed, where a positively correlated bivariate Gaussian distribution denotes the component of reproducible sites and a zero-correlated bivariate Gaussian distribution denotes the component of irreproducible sites.

where $(x,y)$ denotes site ranks of each pair of techniques, $p$ denotes the proportion of the reproducible site component, and $\mu$, $\sigma^2$, and $\rho$ denote the parameters for the positively correlated bivariate Gaussian distribution.

A two-component semiparametric copula mixture model is employed to assess the dependence structure between the ranks of paired techniques. The associations among the variables are modelled by a parametric joint distribution, while the marginal distributions are estimated nonparametrically based on their ranks. The Gaussian copula family is employed to model dependence structures, treating the observed data as being generated from underlying latent variables $(z_x, z_y)$. The underlying latent variables reflect the biological replicates and are assumed to have the same marginal distribution. Therefore, the joint parametric model for copula with latent variables $(z_x, z_y)$ is:

where $K_i \sim \text {Bernoulli} (p)$, $\mu_0 = 0$, $\mu_1 > 0$, $\rho_0 = 0$, and $0<\rho_1<1$.

The IDR value is the posterior probability that a rank pair $i$ denoted as $(x_i, y_i)$ belongs to the irreproducible component:

where $h_0 \sim \text{Normal} \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \right)$ and $h_1 \sim \text{Normal} \left( \begin{pmatrix} \mu_1 \\ \mu_1 \end{pmatrix}, \begin{pmatrix} \sigma_1^2 & \rho_1 \cdot \sigma_1^2 \\ \rho_1 \cdot \sigma_1^2 & \sigma_1^2 \end{pmatrix} \right)$. $F_x$ and $F_y$ are the unknown marginal distributions of ranks $x$ and $y$ respectively. $G^{-1}$ is are the right-continuous inverses of $G$ which is defined as:

where $F_x^{-1}$ and $F_y^{-1}$ are the right-continuous inverses of $F_x$ and $F_y$.

To fit procedure for this model is implemented with a EM algorithm coupled with “pseudo-likelihood” likelihood approach. For detailed fitting procedure, please refer to Li et al. (2011).

3. Threshold selection

By setting a threshold for the IDR value, the irreproducible discoveries can be controlled. In m6AConquer, this threshold is set to 0.05, which means the overlapped sites with posterior probabilities of belonging to the irreproducible component of less than 5% are considered reproducible.

Technical orthogonal integration

IDR analysis is conducted between orthogonal technique pairs. Here, 10 m6A quantitative profiling techniques are classified into 4 groups based on their sequencing principles:

• Chemical-assisted: Uses specific chemicals to deaminate normal A sites but keep m6A modification sites intact

• Enzyme-assisted: Uses enzymes to capture or cleave at m6A modification sites

• Antibody-assisted: Uses the antibody-immunoprecipitation approach to capture RNA fragments with m6A modifications

• Direct RNA sequencing: Deciphers the m6A modification information by detecting distinct current shift patterns compared to their unmodified counterparts.

Technqiues in different strategy groups are considered orthogonal technique pairs. One exception can be the techniques in antibody-assisted and direct RNA sequencing strategy groups. They are not considered orthogonal because the computational detection method in direct RNA sequencing techniques incorporates m6ACE-seq sequencing results as its training labels. m6ACE-seq is classified in the antibody-assisted strategy group.

In each orthogonal technique pair, overall and condition-specific sample integration methods can be employed. The overall integration method is the sum of m6A counts and total counts at each site from all samples within each technique. As for the condition-specific integration method, summation takes place on samples with common cell lines or tissues.

The IDR analysis-based integration results are presented in Genomic Ranges and CSV format in the m6AConquer database. For each orthogonal technique pair, the results contain the genomic locations and IDR values of all reproducible sites, in addition to their m6A methylation ratios, posterior probabilities, and Benjamini-Hochberg (BH) corrected p-values as modelled by each technique. Furthermore, a merged Genomic Range and CSV files are provided, containing the union of reproducible sites identified in all orthogonal technique pairs, along with the names and numbers of supported orthogonal technique pairs. Users can access and download these results from the Download page.

V. Codes

The codes for processing sequencing data is available on GitHub.

The fitting of beta-binomial mixture model and other candidate models can be implemented through OmixM6A R package, which is also available on GitHub.

VI. References

Morgan M, Obenchain V, Hester J, Pagès H (2024). SummarizedExperiment: SummarizedExperiment container. R package version 1.34.0, https://bioconductor.org/packages/SummarizedExperiment.

Ramos M, Schiffer L, Re A, Azhar R, Basunia A, Rodriguez Cabrera C, Chan T, Chapman P, Davis S, Gomez-Cabrero D, Culhane A, Haibe-Kains B, Hansen K, Kodali H, Louis M, Mer A, Reister M, Morgan M, Carey V, Waldron L (2017). “Software For The Integration Of Multi-Omics Experiments In Bioconductor.” Cancer Research, 77(21), e39-42. doi:10.1158/0008-5472.CAN-17-0344, https://cancerres.aacrjournals.org/content/77/21/e39.

Qunhua Li. James B. Brown. Haiyan Huang. Peter J. Bickel. "Measuring reproducibility of high-throughput experiments." Ann. Appl. Stat. 5 (3) 1752 - 1779, September 2011. https://doi.org/10.1214/11-AOAS466

Konstantin Krismer, Yuchun Guo, David K Gifford, IDR2D identifies reproducible genomic interactions, Nucleic Acids Research, Volume 48, Issue 6, 06 April 2020, Page e31, https://doi.org/10.1093/nar/gkaa030